1 Executive Summary

  • Study groups separate well following dimensionality reduction and clustering analysis

  • ~10000 genes (at FDR<0.05) are differentially expressed between MSC[SCI] and MSC[Naive]

  • MSC[Naive] is similar to MSC[in vitro] while MSC[SCI] is different from MSC[in vitro]

  • Genes differentially expressed between MSC[SCI] and MSC[Naive] have functions (GO terms) mainly related to immune system and immune response. Most relevant KEGG pathways for this contrast are related to immune system diseases

  • Heatmap analysis reveals 7 clusters separating MSC[in vitro] from MSC[SCI] and MSC[Naive]. 700 genes (FDR<1e-6) in cluster 2 explain most of the difference between MSC[SCI] and MSC[naive]

2 Project Description

Hypothesis: The transcriptome of bone marrow derived mesenchymal stem cells (MSCs) transplanted into a spinal cord injury environment is altered as compared to MSCs transplanted into uninjured spinal cord.

Experimental setup & sequencing: mCherry+MSCs were transplanted into injured (75 kdyn contusion injury, 24h post SCI) or uninjured spinal cord using a glas capillary pipette. C57BL/6J female mice were used in all experiments. At 7 days post transplantation the mCherry+MSCs were collected using FACS. FACS of MSCs from injured spinal cord was performed using injured spinal cord as reference (i.e. negative gate setting). FACS of MSCs from uninjured spinal cord was performed using uninjured spinal cord as reference. Total RNA, digested of DNase, was isolated and libraries were prepared using SMARTer Stranded Total RNA-Seq Kit - Pico Input Mammalian (Takara). Sequencing was performed 2x125bp in two lane using the HiSeq2500 system and v4 sequencing chemistry (Illumina Inc.) performed by the SNP&SEQ Technology Platform (Stockholm, Sweden).

Data analysis: Data was analysed using edgeR and limma packages, both available through Bioconductor using R (version 3.4.4, Someone to Lean On). Two additional key packages used were ggplot2 and data.table.

3 Data Overview

Table 1. Data set characteristics.

Characteristic Value
Samples (n): 22
Groups (n): 3
Unique ENSEMBL IDs (n): 46078

4 Data Pre-Processing

4.1 Removing genes with low expression

Background: The expression of a gene must reach a certain threshold for it to be translated into a protein. Translation into a protein is a prerequiste for a gene to have any biological function. Thus, genes with low number of read counts across samples are probably not differentially expressed and should be removed. However, a greater sequencing depth (i.e. a larger library size) will result in a higher read count, thus introducing a bias to the DGE analysis. Therefore, prior to filtering, raw counts are transformed into counts per million (CPM) which accounts for the difference in library size.

Transformation: Raw counts are transformed into CPM by dividing each count with the sum of the column (i.e. the library size) and then multiplying by \(1e6\). Log-CPM is calculated by taking \(log_2(count+0.25)\).

Filtering: Genes with low expression are removed from the count matrix by filtering. A gene is defined as highly expressed if \(CPM>1\) (i.e. \(log(CPM)>0\)) for at least three samples.

Fig 1. Density of log-CPM values pre -and post filtering

Fig 1. Figure reports the density of log-CPM for every sample (by color) pre -and post filtering of genes with low expression. The raw read count matrix is filtered based on log-CPM values. Vertical dashed line represents the cut-off (log-CPM=0, CPM=1). The figure shows a distinct shift of the density from below the threshold (Fig 1A) to above the threshold (Fig 1B). Approximately 1/3 of the genes remain post filtering.

4.2 Adjusting gene expression distributions

Background: The read count is affected by: 1) the gene expression and 2) the sequencing depth. The sequencing depth equals the library size. The library size is defined as the sum of counts for each sample (i.e. column sum). The counts per sample represent the relative abundance of each gene. Highly expressed genes can consume a substantial proportion of the library size thus making the other genes seem underexpressed. Therefore, normalization is conducted in order to ensure that the distribution of the expression is similar for each sample. All samples should have a smiliar range and distribution of expression (log-CPM).

Normalization: Scaling factors are calculated using the trimmed mean of M-values (TMM) algorithm. The algorithm finds a set of scaling factors which minimizes the log-fold change between the samples. Scaling factors >1 downscale the counts while scaling factors <1 scale the counts upwards. The effective library size is then obtained by taking the product of the original library size and the scaling factor for each sample respectively. The effective library size is then used in down-stream calculations.

Fig 2. Distribution of log-CPM values pre -and post normalization

Fig 2. Figure reports the distribution of gene expression (log-CPM) for each sample. Fig 2A reports the distribution prior to normalization while Fig 2B reports the distribution following normalization of library sizes using the TMM algorithm. Boxplots are based on all log-CPM values while points represent a random sample of 1e4 genes (due to processing time issues). The difference in the distribution of log-CPM using original and effective library sizes is minor but adjusted for.

5 Dimensionality Reduction

Process of reducing the number of random variables by obtaining a set of principle variables. The number of dimensions of a dataset is reduced without compromising the information content. Two main approaches: 1) feature selection (find a subset of the original variables) and 2) feature extraction (transformation from high-dimensional to low-dimensional space using linear or non-linear techniques).

Fig 3. Variance explained by principal components based on the 500 genes with highest variance

Fig 3. Figure reports the proportion variance explained by each principal component. Fig 3A reports the proportional variance explained by each component while Fig 3B reports the cumulative variance explained by the components. It is obvious that the first principal component explains the absolute majority of the variance while the remaining components explain only a small portion of the variance.

Table 2. Upper and lower bounds (bootstrapped 95 % confidence intervals) for the proportion of variance explained by principal component 1 to 10

PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
Upper bound: 0.958 0.006 0.004 0.002 0.002 0.002 0.002 0.002 0.002 0.002
Lower bound: 0.967 0.008 0.006 0.003 0.003 0.002 0.002 0.002 0.002 0.002

Fig 5. Dimensionality reduction using PCA and t-SNE of the 500 genes with highest variance

Fig 5. Figure reports the samples in low dimensional space following dimensionality reduction with PCA and t-SNE for the 500 genes with the highest variance. Ellipses are added to the samples for easier recognition of the study groups.

Fig 6. Dimensionality reduction using PCA and t-SNE of the 500 genes with highest variance

Fig 6. Figure reports the samples in three-dimensional space following dimensionality reduction with PCA and t-SNE. The third dimension is added in order to visualize the separation when the influential first component/variable is not used. Two-dimensional representations are added on the walls in addition to the three-dimensional representation.

Fig 7. 1000 bootstrap replicates of PCA of the 500 genes with highest variance

Fig 7. Figure reports the three-dimensional representation of 1000 bootstrap replicates of PCA. Left figure reports with three-dimensional coordinates while right figure does not.

Fig 8. Dimensionality reduction using other common methods for the 500 genes with highest variance

Fig 8. Figure 8A-F reports common dimensionality reduction algorithms applied on the dataset.

6 Clustering Algorithms

Fig 9. Hierarchical clustering of samples using 500 most variable genes

Fig 9. Figure reports hierarchical clustering based on the 500 most variable genes. Fig 9A reports the clustering using a dendrogram while Fig 9B reports the same clustering using a circular packing plot.

Fig 10. Clustering of samples using common methods following dimensionality reduction with PCA and t-SNE

Fig 10. Figure reports result of common clustering algorithms implemented on samples following dimensionality reduction using PCA and t-SNE. In the case train-test split was required a 60:40 ratio (13:9 samples) was used instead of a traditional 80:20 ratio due to low number of samples in the test set.

Fig 11. Clustering of 200 random bootstrap replicates for each sample using affinity propagation and K-means clustering

Fig 11. Figure reports three-dimensional representations of clustering using affinity propagation and K-means clustering of 200 random bootstrap replicates of PCA per sample.

Fig 12. Top ten loadings (absolute value) for the first four principal components

Fig 12. Figure reports the top ten loadings (genes) using absolute values for the first four principal components. No label indicates NA value (no symbol available).

Fig 13. Top 10 positive and negative loadings for the first -and second principal component

Fig 13. Figure reports the top 10 positive and negative loadings for the first -and second principal component.

7 Negative Binominal Dispersions by Weighted Likelihood Empirical Bayes

Fig 14. Cox-Reid profile-adjusted likelihood estimated tagwise, common and trended dispersions

Fig 14. Figure reports tagwise, common and trended dispersions.

8 Linear Modelling and Empirical Bayes Moderation Using Precision Weights

Background: The variance of raw counts is not independent of its mean. Furthermore, the variance of raw counts increases with the count size (opposite is true for log-counts).

Voom transformation: Acronym for mean-variance modelling at the observational level. The mean-variance relationship is estimated in the data which is then used for computing precision weights for each gene. The precision weights are then implemented in the linear modelling in order to adjust for heteroscedasticity. Estimation of precision weights is done as follows: 1) log-CPM=\(log_2(count + 0.5)\), 2) a linear model is fitted to each gene, 3) lowess function is fitted to the scatterplot between average log-CPM and sqrt(st.dev), 4) the variance of the log-CPM for each gene is predicted using the lowess trend line 5) the inverse of the variance for each log-CPM value is the precision weight.

Empirical Bayes moderation of standard error: Ranks genes in order of evidence of differential expression.

Fig 15. Mean-variance relationship pre -and post voom transformation

Fig 15. Figure reports the mean-variance relationship pre -and post application of the voom function. Fig 15A reports the average log-CPM against the quarter root of the variance. Fig 15B reports average log-CPM against the \(log_2(st.dev)\). Blue line reports the average \(log_2(st.dev)\). The red line is a linear trend fitted to the black dots. Each black dot represents a gene. Fig 15A illustrates that the variance is decresing when the average expression is increasing. In Fig 15B the dependency is removed and the mean variance is unchanged when the average expression increases.

Fig 16. Hierarchical clustering and circular packing plot of 1000 genes with highest F-value

Fig 16. Figure reports hierarchical clustering of samples based on the 1000 genes with highest F-values. Fig 16A reports a dendrogram while Fig 16B reports a circular packing plot.

9 Differential Gene Expression

Determining significant gene differential expression for a contrast: A simple Bayesian model is used for moderating standard errors across genes (squeeze them towards a common value) and producing moderated t-statistics. These moderated t-statistics are used for significance analysis. Moderated t-statistics have a higher degree of freedom in comparison to usual t-statistics due to the increase in reliability resulting from the smoothening of standard errors. P-values are adjusted for multiple testing using Benjamini and Hochberg’s method to control for false discovery rate (FDR).

Fig 17. Number of differentially expressed genes for each contrast (FDR<0.05)

Fig 17. Figure reports a venn diagram containing the number of differentially expressed genes for each contrast.

Table 3. Number of differentially over -and under-expressed genes for each contrast (FDR<0.05)

MSC[SCI]-MSC[naive] MSC[SCI]-MSC[invitro] MSC[naive]-MSC[invitro]
Downregulated: 5248 5301 909
No change 5148 5296 13711
Upregulated: 4888 4687 664
Sum: 15284 15284 15284

Fig 18. Mean difference -and volcano plot (FDR<1e-6)

Fig 18. Figure 18A reports a mean-difference plot which illustrates the number of over -and under expressed genes. Threshold is set at \(log_2(fold change)\) +/-1 (blue lines). Blue dots represents genes above or below the log-fold change thresholds while red dots represent those genes which are above/below the thresholds and are significantly (p<1e-6) differentially expressed. Fig 18B is a volcano plot which reports the number of significantly (p<1e-6) over -and underexpressed genes (marked with red). Blue dots represent genes which have logFC <-1 or >1 but are not significantly expressed. Data in both plots is for the contrast “MSC[SCI] vs MSC[naive]”.

Table 4. 10 most significantly up -and downregulated differentially expressed genes

Gene log2(fold change) P-value (adjusted) Gene log2(fold change) P-value (adjusted)
Lcp1 -13.37 1.1e-21 Aars 1.61 1.3e-14
Hpgds -9.92 1.1e-21 Hspa9 1.67 1.1e-13
B430306N03Rik -9.17 1.4e-21 Slc7a1 2.09 1.5e-13
St8sia4 -9.10 1.4e-21 Slc7a5 2.52 1.7e-13
Cd48 -9.34 1.8e-21 Txnrd1 1.54 1.9e-13
Cybb -12.95 1.0e-20 Gars 1.48 1.0e-12
Ms4a7 -11.60 1.0e-20 Hells 2.09 1.1e-12
March1 -10.49 1.0e-20 Kitl 2.51 1.3e-12
Ptafr -9.85 1.0e-20 Utp20 1.67 1.2e-12
Ccr1 -10.17 1.2e-20 Urb2 1.47 1.8e-12

10 Gene Ontology and KEGG Enrichment Analysis

Gene Ontology (GO): A major bioinformatics initiative which offers a computational representation of the biological function of genes at molecular, cellular and tissue level. This tool enables one to annotate genes with their function.

KEGG (Kyoto Encyclopedia of Genes and Genomes): KEGG pathway are manually drawn pathway maps which represent the current knowledge within metabolism, cellular processes and many more.

Table 5. GO terms and KEGG pathways (FDR<1e-6)

Type Term Ont N Up Down P.Up P.Down
GO immune system process BP 1566 107 542 1.0e+00 8.2e-116
GO immune response BP 828 50 358 1.0e+00 4.1e-105
GO defense response BP 858 51 348 1.0e+00 2.6e-92
GO regulation of immune system process BP 857 48 328 1.0e+00 9.3e-79
GO positive regulation of immune system process BP 585 30 252 1.0e+00 5.0e-72
GO regulation of immune response BP 454 23 212 1.0e+00 6.7e-68
GO leukocyte activation BP 570 34 238 9.9e-01 1.1e-64
GO cell activation BP 657 39 259 1.0e+00 2.0e-64
GO inflammatory response BP 428 22 193 1.0e+00 1.2e-58
GO response to external stimulus BP 1260 96 376 9.3e-01 9.0e-57
Type Pathway N Up Down P.Up P.Down
KEGG Cytokine-cytokine receptor interaction 149 7 75 9.8e-01 9.1e-21
KEGG Staphylococcus aureus infection 38 0 32 1.0e+00 3.0e-19
KEGG Leishmaniasis 60 0 40 1.0e+00 2.5e-17
KEGG Osteoclast differentiation 114 3 58 1.0e+00 1.5e-16
KEGG Lysosome 112 0 55 1.0e+00 6.7e-15
KEGG Inflammatory bowel disease (IBD) 44 1 31 9.8e-01 9.7e-15
KEGG Hematopoietic cell lineage 56 5 34 5.7e-01 3.9e-13
KEGG Antigen processing and presentation 60 2 35 9.8e-01 9.2e-13
KEGG Tuberculosis 144 7 59 9.8e-01 1.3e-11
KEGG Phagosome 132 3 55 1.0e+00 3.0e-11

11 Unsupervised Clustering 2

11.1 Agglomerative hierarchical clustering with heatmap

Background: Purpose of hierarchical clustering in combination with heatmap is to find subset of genes which explain the difference between study groups.

Hierarchical clustering: An unsupervised statistical learning method which aims at clustering the data in a not pre-determined amount of clusters. In agglomerative hierarchical clustering the tree is built from the terminal nodes (leaves) towards the root. A dendrogram is utilized to displayed the clusterings. A heatmap is added to the hierarchical clustering to enable the identification of gene clusters which might explain the hierarchical clustering.

Fig 19. Hierarchical clustering of samples together with heatmap of significantly differentially expressed genes

Fig 19. Figure reports a heatmap with hierarchical clustering (indicated with dendrograms) using log-CPM values. Only significantly differentially expressed genes are included (and genes with NA symbols were removed).

Table 6. 10 most significantly up -and downregulated differentially expressed genes in cluster two

Gene log2(fold change) P-value (adjusted) Gene log2(fold change) P-value (adjusted)
Lcp1 -13.37 1.1e-21 Gsta1 6.01 2.3e-13
C1qb -12.93 2.9e-19 Cyp1a1 5.74 7.0e-12
Slamf7 -12.01 3.6e-20 Usp17la 5.03 6.1e-10
Ms4a6c -11.86 4.6e-22 Ptchd4 4.92 1.1e-11
H2-Aa -11.63 1.7e-18 Nqo1 4.75 4.1e-11
Xist -11.48 1.3e-17 Gsta4 4.29 5.9e-12
Cd84 -11.22 5.6e-21 Cyp1b1 3.81 3.7e-09
Fcer1g -11.17 8.1e-18 Sema3a 3.58 1.1e-08
Dock2 -11.07 5.1e-17 Gsta2 3.52 5.7e-09
Fgl2 -10.92 1.9e-15 Erfe 3.25 8.9e-10

Fig 20. Chord diagram representing relation between study groups and clusters in heatmap

Fig 20. Figure reports a chord diagram which relates the study groups with the associated clusters.

11.2 Multiscale bootstrap resampling of agglomerative hierarchical clustering

Fig 21. Multiscale bootstrap resampling of hierarchical clustering

Fig 21. Figure reports hierarchical clustering using multiscale boostrap resampling with 1000 replicates. Numbers with red color are the approximately unbiased (AU) p-values while blue numbers represent bootstrap probability (BP). AU >95 and BP >70 are considered to be significant clusters. Red boxes indicate clusters with AU >95.

12 Gene Set Testing

12.1 Molecular Signatures Database

Camera: Performs a competitive gene set test which accounts for inter-gene correlations. Camera evaluates if a set of genes is highly ranked relative to other genes in terms of differential expression.

Barcode plot: This function plots two sets of genes in a ranked list of statistics. Statistics are ranked left to right from smallest to largest. Shaded region in the middle of the plot represents the ranked statistics while the vertical bars reports the positions of the specified subsets. The enrichment worm (line) shows the relative enrichment of the vertical bars in each part of the plot.

13 Summary

1. Annotation: Genes are annotated with gene name using their respective ENSEMBL ID.
2. Transformation: Read count matrix is transformed into log-CPM using original library sizes.
3. Filtering: Read count matrix is filtered using log-CPM values (>0 for at least 3 samples).
4. Normalization: Effective library sizes are calculated using the library sizes for the filtered read count matrix and the trimmed mean of M values (TMM) approach.
5. Transformation: Filtered read count matrix is transformed into log-CPM matrix.
6. PCA: Conducted for the 500 genes with highest variance. Proportional variance explained, MDS and loading plots are created.
7. Design matrix: A dummy matrix which indicates which group each sample belongs.
8. Contrast matrix: Contrasts are the group comparisons of interest.
9. Voom transformation: Estimate precision weights for linear modelling to remove dependency between the variance and trhe mean.
10. Linear modelling: Linear modelling using precision weights followed by an empirical Bayes moderation.
11. Differentially expressed genes: Moderated t-statistics are used for determining significantly expressed genes for each contrast. Results are displayed with venn diagrams, mean-difference -and volcano plot and a summary table.
12. Analysis/interpretation: Using hierarchical clustering, heatmap, gene ontology and KEGG enrichment analysis and gene set analysis the difference between the study groups is sought for.

14 Bibliography

[1] R Core Team (2016). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL htts://www.R-project.

[2] Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., Shi, W., and Smyth, G.K. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43(7), e47.

[3] Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139-140

[4] Law CW, Alhamdoosh M, Su S et al. RNA-seq analysis is easy as 1-2-3 with limma, Glimma and edgeR [version 1; referees: 3 approved]. F1000Research 2016, 5:1408.

15 Setup

This analysis was conducted on:

sessionInfo()
## R version 3.4.4 (2018-03-15)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 16.04.3 LTS
## 
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.6.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.6.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=sv_SE.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=sv_SE.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=sv_SE.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=sv_SE.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
##  [1] grid      parallel  stats4    stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2                            
##  [2] circlize_0.4.3                          
##  [3] ape_5.0                                 
##  [4] igraph_1.2.1                            
##  [5] ggraph_1.0.0.9999                       
##  [6] e1071_1.6-8                             
##  [7] apcluster_1.4.5                         
##  [8] MeanShift_1.1-1                         
##  [9] wavethresh_4.6.8                        
## [10] EMCluster_0.2-10                        
## [11] Matrix_1.2-11                           
## [12] kernlab_0.9-25                          
## [13] mda_0.4-10                              
## [14] class_7.3-14                            
## [15] MASS_7.3-47                             
## [16] Rtsne_0.13                              
## [17] caret_6.0-78                            
## [18] lattice_0.20-35                         
## [19] xml2_1.2.0                              
## [20] XML_3.98-1.3                            
## [21] RCurl_1.95-4.10                         
## [22] bitops_1.0-6                            
## [23] pander_0.6.1                            
## [24] knitr_1.20                              
## [25] pvclust_2.0-0                           
## [26] boot_1.3-20                             
## [27] rafalib_1.0.0                           
## [28] ggrepel_0.7.0                           
## [29] plot3D_1.1.1                            
## [30] ellipse_0.4.1                           
## [31] VennDiagram_1.6.19                      
## [32] futile.logger_1.4.3                     
## [33] gplots_3.0.1                            
## [34] RColorBrewer_1.1-2                      
## [35] colorspace_1.3-2                        
## [36] gridExtra_2.3                           
## [37] cowplot_0.9.2                           
## [38] ggplot2_2.2.1.9000                      
## [39] data.table_1.10.4-3                     
## [40] Mus.musculus_1.3.1                      
## [41] TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.0
## [42] org.Mm.eg.db_3.4.1                      
## [43] GO.db_3.4.1                             
## [44] OrganismDbi_1.18.1                      
## [45] GenomicFeatures_1.28.5                  
## [46] GenomicRanges_1.28.6                    
## [47] GenomeInfoDb_1.12.3                     
## [48] AnnotationDbi_1.38.2                    
## [49] IRanges_2.10.5                          
## [50] S4Vectors_0.14.7                        
## [51] Biobase_2.36.2                          
## [52] BiocGenerics_0.22.1                     
## [53] edgeR_3.18.1                            
## [54] limma_3.32.10                           
## 
## loaded via a namespace (and not attached):
##   [1] backports_1.1.2            plyr_1.8.4                
##   [3] lazyeval_0.2.1             splines_3.4.4             
##   [5] BiocParallel_1.10.1        digest_0.6.15             
##   [7] foreach_1.4.4              BiocInstaller_1.26.1      
##   [9] htmltools_0.3.6            viridis_0.5.0             
##  [11] gdata_2.18.0               magrittr_1.5              
##  [13] memoise_1.1.0              sfsmisc_1.1-2             
##  [15] recipes_0.1.2              Biostrings_2.44.2         
##  [17] gower_0.1.2                dimRed_0.1.0              
##  [19] matrixStats_0.53.1         blob_1.1.0                
##  [21] dplyr_0.7.4                graph_1.54.0              
##  [23] bindr_0.1.1                survival_2.41-3           
##  [25] iterators_1.0.9            glue_1.2.0                
##  [27] DRR_0.0.3                  gtable_0.2.0              
##  [29] ipred_0.9-6                zlibbioc_1.22.0           
##  [31] XVector_0.16.0             DelayedArray_0.2.7        
##  [33] ddalpha_1.3.1.1            shape_1.4.4               
##  [35] DEoptimR_1.0-8             scales_0.5.0.9000         
##  [37] futile.options_1.0.0       DBI_0.8                   
##  [39] Rcpp_0.12.16               viridisLite_0.3.0         
##  [41] units_0.5-1                foreign_0.8-69            
##  [43] bit_1.1-12                 lava_1.6                  
##  [45] prodlim_1.6.1              pkgconfig_2.0.1           
##  [47] nnet_7.3-12                locfit_1.5-9.1            
##  [49] labeling_0.3               tidyselect_0.2.4          
##  [51] rlang_0.2.0.9000           reshape2_1.4.3            
##  [53] munsell_0.4.3              tools_3.4.4               
##  [55] RSQLite_2.0                broom_0.4.3               
##  [57] evaluate_0.10.1            stringr_1.3.0             
##  [59] yaml_2.1.18                ModelMetrics_1.1.0        
##  [61] bit64_0.9-7                tidygraph_1.0.0.9999      
##  [63] robustbase_0.92-8          caTools_1.17.1            
##  [65] purrr_0.2.4                RBGL_1.52.0               
##  [67] nlme_3.1-131               RcppRoll_0.2.2            
##  [69] biomaRt_2.32.1             compiler_3.4.4            
##  [71] tweenr_0.1.5               tibble_1.4.2              
##  [73] stringi_1.1.7              highr_0.6                 
##  [75] psych_1.7.8                pillar_1.2.1              
##  [77] GlobalOptions_0.0.13       rtracklayer_1.36.6        
##  [79] R6_2.2.2                   KernSmooth_2.23-15        
##  [81] codetools_0.2-15           lambda.r_1.2              
##  [83] gtools_3.5.0               assertthat_0.2.0          
##  [85] CVST_0.2-1                 SummarizedExperiment_1.6.5
##  [87] rprojroot_1.3-2            withr_2.1.2               
##  [89] GenomicAlignments_1.12.2   Rsamtools_1.28.0          
##  [91] mnormt_1.5-5               GenomeInfoDbData_0.99.0   
##  [93] udunits2_0.13              rpart_4.1-11              
##  [95] timeDate_3043.102          tidyr_0.8.0               
##  [97] rmarkdown_1.9              misc3d_0.8-4              
##  [99] ggforce_0.1.1              lubridate_1.7.3